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Abstract 


The  Beta-Gamma  transformation  is  described  and  is  used  to  define  a  very 
simple  first-order  autoregressive  Beta-Gamma  process,  BGAR(1).  Maximum 
likelihood  estimation  is  discussed  for  this  model,  as  well  as  moment 
estimators.  The  first-order  structure  is  extended  to  include  moving  average 
processes  and  mixed  first-order  autoregressive,  pth-order  moving  average 
processes.  It  is  shown  that  these  Gamma  processes  are  time-reversible  and, 
therefore,  too  narrow  for  general  physical  modelling.  A  dual  process  to  the 
BGAR(1)  process,  DBGAR(i),  is  introduced,  as  well  as  an  iterated  process 
which  combines  the  Beta-Gamma  process  and  the  GAR(  1 )  process  of  Gaver  and 
Lewis(1980).  Some  properties  of  these  extended  autoregressive  processes  are 
derived.  Several  highly  nonlinear  extensions  of  these  processes  which 
produce  negative  correlation  are  given.  Use  of  the  processes  to  model  a 
sequence  of  times  between  failures  of  a  computer  system  is  described. 

0.  INTRODUCTION 


The  Gamma  distribution  is  used  to  model  a  wide  variety  of  positive 
valued  random  quantities  in  fields  such  as  operations  analysis,  reliability 
theory,  hydrology  and  meteorology.  Thus,  service  time  distributions  and 
interarrival  times  in  queues  are  often  modelled  as  having  Gamma  distribu¬ 
tions,  as  are  wind  velocities  (Hugus,  1982;  Brown,  Katz  and  Murphy,  198^4) 
measured  at  successive  discrete  time  points  or  river  flows  at  successive 
instants  of  time  (Lawrance  and  Kottegoda,  1977).  In  all  these  cases,  the 
measurements  are  taken  serially  in  time  and  are  apt  to  be  serially 
dependent.  Thus,  development  of  time  series  with  Gamma  distributed  marginal 
distributions  and  various  correlation  structures  is  of  great  importance. 

Gaver  and  Lewis  (1980)  showed  that  the  usual  linear  first-order 
autoregressive  equation,  X  «  pX  ,♦£  ,  would  yield  Gamma  marginal 
distributions  if  the  i.i.d.  sequence  {E^}  was  chosen  with  suitable  marginal 
distribution.  This  Gamma  innovation  distribution  has  a  positive  probability 


i 


1 


.. .  vy  Codts 

Avail  and/or 
lipacM 


□  □ 


of  being  zero,  so  that  the  process  (GAR(1))  generates  sample  paths  which 

exhibit  'runs-down'  (as  seen  in  river  flow  data),  but  which  are  "defective". 

The  defect  lies  in  the  fact  that  when  ■  0,  and  ,  are  proportional 

n  n  n“i 

and  p  can  be  estimated  exactly  in  long  enough  time  series.  Moreover,  the 
probability  of  the  defect  is  higher  for  k  small,  which  is  precisely  where 
the  model  is  needed,  since  the  usual  techniques  of  transforming  to  normality 
are  then  questionable  and  probably  undesirable. 

Bernier  (1970)  introduced  the  GAR(1)  model  in  a  hydrological  context 
and  McKenzie  (1982)  introduced  a  multiplicative  Gamma  process  called  PAR(l) 
-  product  autoregression  of  order  one. 

In  Lewis  (1983)  and  Hugus  (1982),  a  simple  linear,  random  coefficient 
model  called  BGAR(1)  was  introduced.  It  is  based  on  the  Beta-Gamma 
transformation  described  in  Section  1. 

The  purpose  of  this  paper  is  to  develop  the  properties  of  this  BGAR(1) 
model  and  to  extend  the  idea  to  moving  average  and  mixed  autoregressive 
structures.  In  particular,  it  is  shown  that  these  processes,  like  the 
Gaussian  ARMA(p,q)  models,  are  time  reversible  and  therefore  are  very 
particular . 

Several  schemes  for  broadening  the  structure  of  Gamma  time  series  are 
given.  In  particular,  a  technique  of  iteration  produces  a  Gamma  autoregres¬ 
sive  process  with  two  structional  parameters  that  can  model,  for  given 
marginal  distribution  and  serial  correlation,  different  kinds  of  sample  path 
behavior.  Some  nonlinear  schemes  that  produce  negative  serial  correlation 
are  also  introduced. 

It  is  important  to  note  the  multiplicity  of  Gamma  processes  which  can 
be  derived  with  given  first-  and  second-order  structure.  Conseqently,  in 
the  absence  of  a  'natural’  structure  such  as  exists  for  Gaussian  processes, 
our  aim  has  been  to  produce  simple  structures,  i.e.,  linear,  additive, 
random  coefficient  processes. 

Finally,  a  series  of  times  between  failures  of  digital  computers 
(Lewis,  196^4)  is  analyzed  and  fitted  with  the  model.  The  data  is  serially 
correlated  with  a  marginal  distribution  which  is  more  skewed  than  an 
exponential  distribution.  Although  this  data  is  known  to  be  generated  by  a 
branching  Poisson  process  (cluster  process),  the  Gamma  model  is  much  simpler 
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and  much  more  tractable  than  the  cluster  process,  and  provides  an  adequate 
representation  for  most  purposes. 


1 .  PRILXHXNARXIS 


In  what  follows  we  will  use  B(m,n),  or  simply  B  when  the  parameteriza¬ 
tion  is  clear,  to  stand  for  a  Beta  random  variable  with  parameters  m  >  0  and 
n  >  0,  denoted  by  Beta(m,n).  The  probability  density  function  for  a 
Beta(m,n)  random  variable  is 


fg(x;m,n) 


r(m*n) 

r(m)r(n) 


m-1 .n-1 
X  (1-x) 


0SxS1;m>0;n>0, 


(1.1) 


where  r(*)  is  the  complete  Gamma  function. 

We  will  denote  by  {B(m,n),  B’(m' ,n’ ) , • • • }  an  i.i.d.  sequence  of  vector 
random  variables  whose  components  are  independent  Beta  random  variables. 

Let  G(k,6)  stand  for  a  Gamma  random  variable  with  shape  parameter 
k  >  0,  and  rate  parameter  B  >  0,  denoted  by  Gamma(k,6).  The  probability 
density  function  for  a  Gamma(k,6)  random  variable  is 

-k  k-1  -fix 

fG(x;k,6)  -  - ,  X  ^  0;  B  >  0;  k  >  0.  (1.2) 

We  will  denote  by  (G^(k,B)},  an  i.i.d.  sequence  of  Gamma  variates. 

A  Gamma (k, 6)  random  variable  has  moments 

E(G)  -  k/6  -  w;  Var(G)  -  k/6^;  C(G)  -  s.d(G/u)  -  ,  (1.3) 


where  C(G)  is  the  coefficient  of  variation,  and  Laplace-Stieltjes  transform 

Lq(u)  -  E(e"''^)  -  (g-T-;;)''- 

The  Gamma  variable  is  sometimes  parameterized  in  terms  of  the  parameter 
p  •>  E(G).  This  is  useful  in  statistical  work,  since  the  mean  is  a  multi¬ 
plicative  parameter  and  G  can  be  written  as  a  unit-mean  Gamma  variate,  G*, 


times  u,  i.e.,  G  ■  pG*.  However,  in  what  follows,  we  will  use  the  fact  that 
two  independent  Gamma  variates  with  the  same  0 “parameters,  but  possibly 
different  shape  parameters,  add  to  give  another  Gamma  variate 

G"(k^k',6)  -  G(k,B)  ♦  G'(k',a).  (1.5) 


The  result  is  not  true  if  the  Gamma  variates  have  the  same  mean  but 
different  shape  parameters. 

The  Gamma  family  of  random  variables  include  the  Exponential(k«1 ) , 
Erlang(k  integer)  and  Chi“Square(k“r/2,  r-1,2,***;  B-2)  random  variables. 

Gamma  and  Beta  variates  are  intimately  related  and  two  of  their 
properties  will  be  used  throughout  this  paper. 

(i)  A  Beta(m,n)  variate  may  be  generated  as 


B(m,n) 


G»(m.6) _ 

G'(m,8)  ♦  G-(n,6)’ 


(1.6) 


where  G'(m,6)  and  G"(n,6)  are  independent.  Furthermore,  the  ratio  B(m,n)  is 
independent  of  the  denominator  G'(m,6)  +  G"(n,6)  ■  G(m*n,6)  and  this 
property  characterizes  the  Gamma  random  variable  (Johnson  and  Kotz,  1970a). 

(ii)  The  Beta  Gamma  transformation.  Multiplying  a  G(m+n,B)  random 
variable  by  an  independent  B(m,n)  random  variable  gives  a  G(m,6)  random 
variable 


G(m,6)  -  B(m,n)G(m*n,6) .  (1.7) 

Thus  one  can  reduce  the  shape  parameter  of  a  Gamma  random  variable  (multiply 
by  a  Beta  random  variable),  as  well  as  increase  it  (add  an  independent  Gamma 
variate,  as  at  (1.5)).  A  heuristic  argument  for  the  result  (1.7)  is  that  if 
we  wanted  to  perform  the  operation  (1.7)  on  a  computer,  we  could  first 
generate  G'(m,6)  and  G"(m,6)  to  form  the  ratio  (1.6)  to  obtain  the  B(m,n) 
variate.  There  is,  however,  no  need  to  generate  G(m+n,6)  in  (1.7),  we  can 
use  G'(m,6)  ♦  G"(n,6),  which  is  Independent  of  the  Beta  variate. 
Multiplication  then  gives  G(m,6)  ■  G'(m,6). 
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(iii)  A  formal  proof  of  the  Beta-Gamma  transformation  is  a  special 
case  of  the  following  Lemma,  which  will  be  used  in  the  sequel. 

kemma.  Let  X(k,6)  be  a  Gamma  random  variable  and  let  B(kp,kp)  be  a 
Beta  random  variable,  which  is  independent  of  X(k,B),  with  p  -  1-p 
lying  in  (0,1).  Then 


,  -(v>Bu)X,  _  f_6_^kp  (  B  ^kp 
^  '  '■B+v^  '•B*v^u^ 


V  a  0,  u  2  0. 


(1.8) 


When  V  -  0,  this  result  proves  the  Beta-Gamma  transformation  in  the 

form 


X(kp,B)  •  B(kp;kp)X(k,B). 


(1.9) 


Proof .  Using  (1.4)  and  conditioning  on  B,  we  get 


where  we  have  used  the  finiteness  of  the  expectation  to  take  the  expectation 
inside  of  the  binomial  expansion. 

ff  —  -  - 

Now  since  E(B  )  -  B(kp*«,,kp)/B(kp,kp)  -  r(kp+a)  r(kp)r(k )/{  r(k  +  £)  r(kp) r(kp) } 
we  have  that 


fk  +  il-1 
(  i- 


E(B*') 


fk+l-l 

I 


) 


r(kp^£)r(kp)r(k) 

r(k+il)r(kp)r(kp) 


r(k*i!,)r(kp^ii)r(kp)r(k) 

r(k)r(«.^i)r(kp)r(kp)r(k+t) 
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Thus 


I  ("T’x- .  I  jS-)' 


£■0 


£-0 


6+v 


1 

1.  ® 

r"! 

[  B  j 

1,.  u  )  ■ 

B^v  J 

[B+v+u 

1 1 

and  therefore 


E|,-(»*Bu)X, 


/  6  xkpx  6  ^kp 
'■6+v*u'^ 


which  was  to  be  proved. 


2.  THE  FIRST-ORDER  AUTOREGRESSIVE  PROCESS,  BCAR(1) 


2.1.  Construction  of  the  Process. 

Using  the  Beta-Gamma  transformation,  we  can  construct  a  very  simple 
first-order  autogressive  process  {X^((<,B,p)}  with  Gamma(k,6)  marginal 
distribution  and  a  single  parameter,  p,  that  describes  the  dependency 
structure  of  the  process.  We  have 

X^(k,6,p)  -  B^(kp,kp)X^_^ (k,6,p)  ♦  B^(kp,kp)G^(k , B)  (0  <  p  <  1)  (2.1) 

-  B^(kp,kp)X^_^ (k,B,p)  +  Y^(kp,B)  n»0,±1,...,  (2.2) 

where  {y^(kp,B)}  are  i.i.d.  Gamma(kp,B)  variates  independent  of  the 
(B^(kp,kp)}  sequence.  If  X^_^(k,B,p)  has  a  Gamma(k,B)  marginal 
distribution,  then  multiplying  by  B^(kp,kp)  reduces  it  to  a  Gamma(kp,B) 
variate  and  adding  the  innovation  variable  Y^(kp,B)  creates  the  Gamma(k,e) 
variate,  X^(k,6).  The  alternate  form  (2.1)  shows  the  process  as  a 
transformation  of  an  i.i.d.  Gamma(k,B)  sequence,  but  clearly  generation  on  a 
computer  would  be  done  with  (2.2).  In  the  sequel,  we  will  drop  the 
parametric  notation  where  no  confusion  is  possible. 
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It  is  01620"  that  taking  to  be  a  Gamma(k,6)  variate  will  start  the 
process  in  a  stationary  mode.  Also,  the  process  is  Markovian  by 
construction. 

2.2.  Serial  Correlation. 

It  is  easily  established,  using  moments  of  Beta  variables,  that 

p(r)  -  Corr(X  ,X  )  «  r  -  0,±1,±2,**‘  .  (2.3) 

n  n~r 

Thus,  in  the  three  parameter  process,  the  parameters  k  and  6  describe  the 
marginal  distribution  of  the  process  and  p  independently  describes  the 
correlation  structure.  Note  that  since  the  process  is  only  defined  for 
0  S  p  <  1  the  correlations  are  non-negative. 

2.3.  Joint  Laplace-Stieltjes  transform. 

In  the  stationary  process  (X  },  let  L  (u,v)  denote  the  joint 

n  A  •  X  a 

n’  n-1 

Laplace-Stieltjes  transform  of  the  adjacent  variables  X  and  X  ,  .  Then  we 

n  n-1 

have 


Ly  y  (u,v)  -  E[exp{-X^u-X„  , 
X  ,x  .  n  n- 1 

n’  n-1 


V}] 


E[exp{-B^X^  ,u-Y  u-X^_,v}] 
n  n-1  n  n-1 


(2.H) 


where,  in  the  last  step,  we  have  dropped  the  indices  n  and  n-1  because  of 
stationarity  and  have  used  the  assumed  independence  of  and  X^_^  to  write 
the  expectation  as  the  product  of  two  expectations. 

Now  the  second  term  is  evaluated  in  the  Lemma  of  Section  1  and  we  have 
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“•x  .X 

n’  n-l 


f  B  ^kpx  B  xkpr  B  ^kp 
'6'*’V+u'^ 


(  6  .  B  ^kpf  B  ^kp 

'B'*’U  B+v'  ^B+v+u' 


(2.5) 


Since  this  transform  is  symmetric  in  u  and  v,  the  joint  distribution  of 

and  is  symmetric.  Also,  since  the  joint  distribution  of  any  set  of 

X^'s  can  be  obtained  from  (2.5)  and  the  marginal  Gamma  distribution,  the 

Beta-Gamma  process  is  tine-reversible. 

Note,  too,  that  we  have  directly  from  the  defining  equation  (2.2)  that 

the  regression  of  X  on  X  ,  -  x  is  linear: 

n  n- 1 

^^^nl^n-l’^)  “  *  (1-p)k/B  -  px  +  (l-p)iJ.  (2.6) 


The  time-reversibility  of  the  process  shows  that 


E(Xn-^  |X^-y)  -  py  +  (1-p)m. 


(2.7) 


2.^.  Convergence  to  a  Gaussian  AR(i)  Process. 

1  /2 

As  k  gets  large,  the  standardized  Gamma(k,1)  variate  X’  «  (X-k)/k 

converges  to  a  standardized  Gaussian  variate.  To  prove  this,  consider  the 

pair  X^  and  X^_^  in  the  BGAR(1)  process.  The  joint  characteristic  function 

for  the  standardized  variables  X’  and  X'  ,  is,  using  (2.5), 

n  n-i 


0^,  (s,t)  -  E{exp(-isX'-itX’  ^ ) } 

■'n’  *n-l 


Taking  logarithms  and  expanding  in  powers  of  k 


-1  /2 


gives 


(s,t)  -  ik^^^(s*t)  -  kp{isk  ^ ^^-( is)^/(2k ) ♦0(k 

^  '  —  -1 /2  2  -■^/2 

-  kp{ itk  -(it)  /(2k)^0(k  ^  )) 

-  kp{  ik~’^^(s+t)-(i)^(s*t)^/(2k)+0(k~^''^)l 
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and  as  this  converges  to 


(s,t)  -  -  |•^s*+t*♦2pst}.  (2.8) 

^n’‘'n-1 

Thus,  since  the  process  is  Markovian,  the  BGAR(1)  process  is  equivalent  to  a 
Gaussian  AR(1)  process  when  k  becomes  large. 

2.5.  Additivity  and  the  GAR(1)  Process. 

Gaver  and  Lewis  (1980)  showed  that  the  usual  linear  stochastic 
difference  equation 


X 

n 


PX  - 


n-1 


(2.9) 


will  give  a  process  (GAR(1))  with  Gamma(k)  marginal  distributions  if  E  is 

n 

the  Gamma  innovation  variable  GI^(k,p)  with  Lapl ace-St i el t j es  transform 
{ ( 1 +ps)/(  1 +s)  where  0  5  p  <  1.  This  variable  can  be  simulated  by  methods 
given  by  Lawrance  (1982)  and  McKenzie  (1986). 

Note  that  the  GAR(1)  process  is  a  linear  additive  (constant 
coefficient)  process  and  adding  two  independent  GAR(1)  processes  with  the 
same  6  and  p  values,  say  {X*(k^)}  and  {X**(k2)}  gives  a  new  GAR(1)  process 
with  shape  parameter  K^+k^  and  dependency  parameter  p.  This  is  not  true  for 
the  BGAR(1)  process  which  is  a  random  coefficient,  linear  additive  process. 
The  process  obtained  by  addition  is  a  process  with  Gamma  marginals  and 

I  r  I 

correlation  structure  p(r)  -  p'  but  it  is  not  even  Markovian.  Additional 
differences  between  the  processes  are  that  while  the  BGAR(1)  process  is 
time-reversible,  the  GAR(1)  process  exhibits  'runs  down'.  In  fact,  it  is 
degenerate  in  the  sense  that  the  innovation  variable  GI^  is  zero  with 
probability  p*^.  Thus,  we  get  ”  P  with  probability  p  ,  and  p  can  be 
estimated  exactly  in  long  enough  series  (Gaver  and  Lewis,  1980).  This 
degenerate  behavior  is  not  exhibited  by  the  BGAR(1)  process.  A  method  for 
combining  two  processes  to  obtain  a  broader  process  is  described  in  Section 


2.6.  The  Conditional  Density  for  X^,  Given  -  y. 

From  the  definition  (2.2),  we  have  that 


P{X^  i  “  P{B^(kp,kp)z  *■  Y^(kp,6)  ^  x} 

-  P{Y  (kp,6)  ^  X  -  B  (kp,kp)z). 
n  n 


Now,  by  definition,  Y^(kp^,8)  is  independent  of  and  of  B^(kp,kp).  Thus 
conditioning  on  B^  and  differentiating  with  respect  to  x  yields  the 
conditional  density  for  X^,  given  X^_^  -  z,  as 


f„  (x,y) - - 

n'^n-1  r(kp){r(kp)}" 


x  ^  0,  y  >  0 


»  jVP“'(1-w)‘<P-Vx^yw)e'^y”dw, 
0 


(2.10) 


where 


L  * 


1  if  X  ^  y, 

x/y  if  X  <  y. 


(2.11) 


and,  for  simplicity,  the  scale  parameter  6  has  been  set  equal  to  one. 

This  density  is  a  continuous  function  of  x  and  absolutely  continuous 
except  where  x  -  y.  Its  utility  is  in  obtaining  maximum  likelihood 
estimates  of  the  parameters  6  (or  p),  k  and  p.  This  is  discussed  in  the 
next  subsection. 


2.7.  Moment  and  Maximum  Likelihood  Estimates  in  BCAR(1). 

There  are  "natural”  moment  estimators  for  the  three  parameters  in  the 
BGAR(1)  model,  namely  the  mean,  u,  (or  g) ,  the  shape  parameter,  k,  and  the 
first-order  serial  correlation  coefficient,  p.  From  (2.3)  and  (1.3),  these 
estimates  are,  from  a  sample  (observed  time  series)  of  length  n. 
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where  S*  is  the  sample  variance,  and 

n-1 

?5  -  I  (X  -X)(X  -X)/{(n-1)S*}.  (2.1iJ) 

i-1  ^  ^  ' 

n 

The  variance  of  p  is  var(X)[l+2  { 1 -(r/n) )  p""  ] /n  and  nVar(u)  is 

r«1 

asymptotically  equal  to  p*  ( 1 +p )/{ ( 1-p)k} .  More  general  properties  of  the 

estimates  R  and  p,  however,  are  hard  to  derive.  But,  their  distributions 

are  independent  of  the  scale  parameter  p. 

In  Table  2.1,  we  give  the  simulated  standard  deviation  and  bias  of  the 

estimates  R  and  p  for  values  of  k  -  0.25,  0.75,  1.0,  2.0  and  p  ■  0.25, 

0.60,  0.90  for  various  values  of  n.  Note  that  the  values  of  n  differ  with 

1  /2 

p,  since  the  "equivalent  sample  size",  n' -nC  ( 1  ♦p)  /  { ( 1 -p)k  }  ]  ,  is 

different.  Here  n'  is  the  sample  size  which  would  be  needed,  for  given  p, 
to  achieve  the  same  variance  for  p  as  in  the  p«0  (independence)  case.  The 
simulation  study  was  performed  with  the  SUPER-SIMTBED  program  of  Lewis,  et 
al.  (1985). 

Two  other  properties  of  the  process  may  be  useful  in  validating  the 
BGAR(1)  model  from  data. 

The  first  is  that  difference  of  successive  values  in  the  time  series, 

D  ■  X  ^  -X  ,  have  an  l-Laplace  distribution  (Dewald,  1985).  This  result 
n  n*1  n 

comes  from  {2.^), 
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Simulation  study  of  the  standard  deviation  and  bias  of  the  moment  estimator,  k,  of  the  shape  parameter,  k,  in  a 
Beta-Gamma  process  for  various  values  of  the  parameters  k  and  p,  and  sample  sizes  n.  In  each  cell,  the  quantity 
in  parenthesis! s  the  estimated  bias  of  k  and  the  top  quantity  is  the  estimated  standard  deviation  of  k. 
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Simulation  study  of  the  standard  deviation  ana  bias  of  the  moment  estimator,  p,  of  the  correlation  parameter,  p,  in 
Beta-Gamma  process  for  various  values_of  the  parameters  k  and  p,  and  sample  sizes  n.  In  each  cell,  the  quantity  in 
parentheses  is  the  estimated  bias  of  p  and  the  top  quantity  is  the  estimated  standard  deviation  of  p. 


after  converting  to  characteristic  functions  and  setting  v 
the  characteristic  function 


•  -u. 


We  have 


1 


“  4  .X 

n  n  n-i 


This  is  the  characteristic  function  of  an  l-Laplace  random  variable  with 
il  -  kp;  the  distribution  is  symmetric  about  zero.  It  goes  to  a  Normal 
random  variable  as  t  but  for  t  ^  1 .  the  density  function  is  not 
absolutely  continuous  at  zero.  In  fact,  for  t  i  1/2,  the  density  is 
infinite  at  zero.  The  fact  that  the  difference  has  median  value  0, 
irrespective  of  the  value  of  p  or  6,  can  be  useful  in  validating  the  model. 

Consider,  now,  ratios 


4  - 

But  Y  (kF,6)  and  X  are  independent  Gamma  variates,  so  that,  if  k>l 
n*1  n 

E(R^)  -  (kp.kp)  ♦  (kp,B)/X^} 

-  p  +  (kp/6)/{(k-1)/6)  -  1  ♦  p/(k-1). 


Higher  moments  can  also  be  obtained. 

These  results  could  be  of  use  in  validating  the  model.  Another 
possibility  for  validating  models  is  the  higher-order  residual  analysis  of 
Lawrance  and  Lewis  (1986). 


Joint-maximum  likelihood  estimates  for  k  and  p  (and  perhaps  p)  can  be 

obtained  from  the  conditional  density  (2.10)  and  the  formula  for  the  joint 

density  of  X  ,X  X,  is 

■'  n  n-1  1 


f(x„;x  ; 
n  n-1 


’’‘l^  ■“  ^X  lx  , 
n'  n-1 


,|x  ‘VrV?’"  • 

n-1 '  n-2 


■“■x,  <’'!>• 


(2.16) 


where  f„  (x.)  is  a  Gamma (k,1)  density. 

1 

Hugus  (1982)  used  (2.16)  to  obtain  joint-maximum  likelihood  estimates 
for  k  and  p.  Three  cases  are  shown  in  Figure  2.1.  Generally,  when  k  is 
greater  than  1,  moment  estimates  of  k  and  p  are  quite  efficient,  which 
agrees  with  results  for  independent  Gamma(k,1)  variates  (Bartlett  and 
Kendall,  19*16).  However,  when  k  is  less  than  one,  the  maximum  likelihood 
estimates  become  much  more  efficient  than  the  moment  estimators. 

The  moment  estimators,  u,  k,  p,  given  at  (2.12),  (2.13)  and  (2.1*1)  serve 
as  good  starting  points  for  numerical  evaluations  to  find  the  maximum 

AAA 

likelihood  estimates,  u,  k,  p,  of  k  and  p.  Techniques  for  the  numerical 
integration  of  (2.10)  are  given  in  Hugus  (1982). 

3.  MOVING  AVERAGE  AND  HIGHER  ORDER  AUTOREGRESSIVE  STRUCTURES 


The  Beta-Gamma  transformation  can  be  used  to  generate  dependency 
structures  other  than  first-order  autoregressive  structures  for  Gamma 
distributed  time  series.  Several  of  these  structures  are  given  in  this 
section. 


3.1.  The  First-Order  Moving  Average  Process,  BGMA(1). 

The  first-order  moving  average  process,  the  BGMA(1),  is  constructed  in 
essentially  the  same  way  as  the  BGAR(1)  above.  If  {G^}  is  a  sequence  of 

i.i.d.  Gamma(k,6)  random  variables  and  {B^}  is  an  independent  sequence  of 

i.i.d.  Beta(ka,ka)  random  variables,  where  a  -  1-a,  we  define  the  (backward) 
BGMA(I)  process  {X^}  by 
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Figure  2. 1 

Simulations  of  joint  maximum  likelihood  estimates  for  k  and  p  and  joint  moment 
estimates  for  k  and  p,  for  three  different  sets  of  values  of  the  parameters. 
Ten  replications  for  each  case.  Top  left  figure:  p-0.75,  k-ii.O;  top  right 
figure:  p-0.25,  k-i4.0;  bottom  figure:  p-0.75,  k-0.75.  The  symbol  (♦)  refers 
to  estimates  of  p;  the  symbol  (o)  refers  to  estimates  of  k.  Large  symbols  are 
joint  estimates;  small  symbols  are  marginal  estimates. 


16 


0  S  a  ^  1. 


(3.1) 


G  +80 


n  n-1 


-1 » 


Evidently,  {X^l  is  a  stationary  process  and  and  are  independent 

if  |r|  >1.  The  marginal  distribution  of  X^  may  be  derived  by  noting  that 

the  right-hand  side  of  (3.1)  is  the  sum  of  a  G(k,6)  random  variable  and  an 

independent  G(l<a,6)  random  variable.  Thus,  X^  is  a  G{k(1^o),6}  random 

variable.  This  process  has  the  same  structure  as  the  usual  Gaussian  MA(1) 

process,  except  that  here  the  coefficient,  B^,  is  a  random  variable  rather 

than  a  constant.  An  immediate  effect  of  this  construction  is  that  the 

observed  and  innovation  processes,  {X  }  and  (G  },  respectively,  have 

n  n 

different  Gamma  marginal  distributions.  This  is  in  contrast  to  the 
structure  of  the  EARMA  processes  (Lawrance  and  Lewis,  1980),  where  it  was 
deliberately  arranged  that  they  should  have  the  same  distribution.  However, 
as  we  shall  see  shortly,  this  disparity  in  marginal  structure  has  some 

V 

advantages. 

From  the  viewpoint  of  modelling,  it  is  more  useful  to  determine  the 

parameters  of  the  innovation  process  in  terms  of  those  of  the  observed 

process.  For  this  reason,  we  reparameterize  (3.1)  slightly.  We  consider 

{G  }  to  be  i.i.d.  Gamma { k / ( 1 +o) , 6 )  r.v.s.  and  {B  }  to  be  i.i.d. 
n  _  n 

Beta{ka/( 1 *0) ,  ka/(1+o)}  random  variables.  This  yields  an  observed  process 

(X  }  which  is  Gamma(k,6). 
n 

We  may  note  that  if  we  write  p  •  Pj((1)  ■  o/(1+a),  then  G^  is 
Gamma(kp,6),  the  same  innovation  process  as  for  the  BGAR(1)  process. 


3.2.  Autocorrelation  Function  for  the  BGMA(l)  Process. 

The  autocorrelation  function  for  the  moving  average  process  may  be 

determined  directly.  Thus,  Cov(X  ,X„  ,)  -  Cov(B  G  ,  ,G  ,)  -  E(B)Var(G), 

n  n-1  n  n-i  n-i 

and  so 


p,(r) 


g 

1>a’ 

0, 


|r|  -  1, 
|r|  >  1. 


(3.2) 


For  |k|  ■  1,  the  attainable  range  of  correlation  is  0  S  p^(1)  i  0.5,  which 
is  the  full  possible  range  of  positive  correlation.  This  is  because,  for  a 
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first-order  moving  average,  |Pjj(1)|  ^  0.5;  see,  e.g.  Hugus  (1982).  The  fact 
that  the  whole  positive  range  is  available  is  important,  because  it  is  in 
contrast  to  the  EMA  models  (Lawrance  and  Lewis,  1977;  1980),  where 
correlation  is  bounded  above  by  0.25.  The  greater  flexibility  in  (3.2)  is  a 
result  of  the  innovation  and  observed  processes  having  different 
distributions.  Since  for  k  -  1,  the  Gamma  distribution  is  an  exponential 
distribution,  the  BGMA(1)  process  is  then  a  broader  first-order  exponential 
moving  average  process  than  the  EMA(1)  process. 


3.3.  Joint  Distributions. 

The  bivariate  Laplace  transform  of  (X  ,X  )  can  be  derived  by  using 

I  n 

(1.8).  Thus,  again,  using  the  notation  a-1-a,  we  have 

L(u,v)  -  E{exp(-uX  -vX  )  -  E{exp(-uG  -uB  G  -vG  -vB  G  ,) 
n+1  n  n+l  n+l  n  n  n  n-l 

-  Lq(u)E{ Lq(v*uB) )Lg^(v) 


[6+uj  [6+vJ  [6+u+vJ 


tka 


ka 


(3.3) 


using  the  Lemma  above.  This  has  exactly  the  same  form  as  the  joint 

transform,  (2.5),  of  BGAR(1)  process  with  o  corresponding 

to  p.  This,  too,  corresponds  to  the  behavior  in  Gaussian  processes,  where 

the  joint  distributions  for  the  autoregressive  and  the  moving  average 

processes  have  the  same  form  and  differ  only  in  their  autocorrelation 

functions.  An  immediate  consequence  is  that  the  conditional  distribution  of 

X  given  X  ,  and  X  given  X  ,  ,  have  exactly  the  same  form  as  for  the 
n* I  n  n  n+i 

BGAR(I).  We  note  the  somewhat  unusual  result  for  a  non-Gaussian  process 

that  regression  is  linear  in  both  directions,  even  though  the  process  is  a 

moving  average.  In  fact,  E(X  |x  -x)  -  E(X  |X  -x)  -  ax^ka/6  -  ax^oiE(X). 

n*i  '  n  n '  n+i 

The  joint  Laplace  transform  of  any  finite  set  of  consecutive 

observations  can  be  obtained  by  the  procedure  that  yielded  (3.3).  Thus,  the 

joint  transform  of  (X  ,X  ,,•••, X  ^, )  is  given  by 

n  n-l  n— r*! 
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L(u^ .u^. 


ka 


.Ur) 


■ika  r-1 

j  *  n 

i-2 


ko  r 

X  n 

i-2 


(3.11) 


Note  that  this  is  not  the  r-dimensional  transform  for  the  BGAR(1)  process. 
Equality  holds  for  only  r  -  2. 

One  consequence  of  (3.^)  is  that,  since  L( u^  .u^  ,  •  •  •  , u^ )  - 
L(u^ ,u^_^ , •  •  • , u^ ) ,  the  process  is  time-reversible. 

The  bivariate  Gamma  distribution  whose  transform  is  given  at  (3.3)  is 
well  known  (Johnson  and  Kotz,  1970b,  p.  219)  and  is  called  by  Ghirtis  (1967) 
the  double  Gamma  distribution.  Since  the  multivariate  form  of  this 
bivariate  Gamma  distribution  arises  as  the  individual  sums  of  m  independent 
Gamma  variates  with  a  common,  independent  Gamma  variate,  it  is  doubtful  that 
triples,  say  X  X  .  ,  X  ,  in  the  BGAR(1)  process  would  have  this 

multivariate  distribution.  In  fact,  (3.*^)  shows  that  this  is  not  so  for  the 
moving  average  process. 

Another  result  that  we  can  immediately  derive  from  the  joint  transform 

is  the  distribution  of  the  sum  of  n  consecutive  observations.  This  has  a 

r 

particularly  simple  form  for  the  BGMA(1)  process.  If  T  -  EX.  then 

"  i-1  ^ 

L.j.(u)  -  L(u  ,u,  •  •  •  ,u) ,  which,  from  (3.9),  is  given  by 


L.j.(u) 


k {r-2(r-1 )o) r  ^  3k(r-1)a 
[B*2uj 


(3.6) 


Further,  since  0  ^  a  S  0.5,  we  can  rewrite  (3.5)  in  terms  of  random 
variables  as 


T^  -  G[k {r-2( r-1 )a) , 6]  +  2G{k (r-1  )a, 6} , 
where  the  two  Gamma  random  variables  are  independent. 

3.9.  Higher  Order  Moving  Average  Processes. 

Higher  order  moving  average  processes  may  be  constructed  by  extending 
the  BGMA(1)  in  an  obvious  way.  Thus,  the  GBMA(q)  is  given  by 
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B  .G  . 
n,i  n-i 


(3.6) 


X  -  G 
n  n 


where  {G  )  is  a  sequence  of  i.i.d.  Gamma{k/(1+  Za.),^)  random  variables  and 
n  ^ 


{B  .  }  is  an  independent  sequence  of  i.i.d.  vector  random  variables  with 

q  q 


3  ..  i  •  1.2,***,q,  being  a  {ka . /(I ♦Za. ) ,  ka./(1*Za.)}  random  variable.  In 


this  case,  {X^}  is  a  stationary  Gamma(k,6)  process  with  X^,  independent 


when  |r|  >  q.  The  autocorrelation  function  is  given  by 


p/r) 


M  '  ;  I  M  , 

a  *  I®-®..  )/fl+  la. 

"  j-' i  j-1  J 


r-1 ,2, . . .  ,q 


r  >  q, 


which,  again,  is  a  close  analogue  of  the  usual  autocorrelation  function  for 
the  Gaussian  MA(q)  process.  The  major  difference  is  that  all  the 


correlations  are  non-negative. 


H.  THE  MIXED  AUTOREGRESSIVE,  MOVING  AVERAGE  PROCESS,  BGARMA(1,1) 


A  more  complicated  dependence  structure  in  Gamma  distributed  variables 
that  is  the  analog  of  the  usual  linear  ARMA(p,q)-type  process  is  now  given. 


^4.1.  Structure  of  the  Mixed  Process.  BGARMA(1,1). 


We  can  construct  an  ARMA-type  process  with  a  Gamma  marginal  distribution 
by  combining  the  two  first-order  processes  we  have  discussed  above.  For 
convenience,  we  write  each  in  a  slightly  different  form.  The  moving  average 


component  is  given  by 


(  -  Y  B  G  , 

n  n-1  n  n 


where  {G  },  {B  }  are  as  given  in  Section  3.1  above,  i.e.,  independent 
n  n 


sequences  of  i.i.d.  Gamma { k / ( 1 *0 ) , 6 }  random  variabes  and  i.i.d. 


Beta{ka/(  1  ■•■a) ,  ka/(1+a)}  random  variables,  respectively.  Notice  that  is 

the  coefficient  of  G  in  (^.1),  whereas  it  was  associated  with  G  ,  in 

n  n-1 

(3.1).  Clearly,  this  change  will  make  no  distributional  difference  and,  as 
we  shall  see,  renders  the  parameters  of  the  ARMA  model  more  readily 
interpretable.  The  sequence  lY^l  is  generated  from  a  BGAR(1)  process  given 
by 


Y  -AY  +  A’G  . 
n  n  n-1  n  n 


('4.2) 


The  process  {G  }  is  as  above  and  {A^}  and  {A^l  are  independent  sequences  of 
i.i.d.  Beta{kp/(  1  ♦&)  ,kp/(1  ^a) }  and  i.i.d.  Beta  I k p/ ( 1 +a) ,kp/ ( 1 +a) }  random 
variables,  respectively.  If  Y^_^  is  G{  (k/(1 +a)  ,6) ,  then  so  also  is  Y^  and 
the  required  stationary  process  results.  The  product  A^G^  is,  of 

course,  the  Gamma{kp/ ( 1 +a) ,6}  random  variable  used  as  the  innovation  process 
in  Section  2.1  above,  but  is  written  in  this  form  here  to  make  explicit  the 
dependence  on  the  innovation  sequence  G^. 

44.2.  The  Autocorrelation  Function  of  the  BGARMA(1,1)  Process. 

The  autocorrelation  function  of  {X^}  may  be  derived  in  the  usual  way. 

Thus,  Cov(X^,X^  )  -  Cov(Y^_, ,Y„  .)  *  Cov(Y^_, ,B^  G  ). 

n  n-r  n-1  n-r-1  n-1  n-r  n-r 

Now,  Cov(Y  „  ,)  -  p''var(Y),  and  Cov(Y  ,B  G  )  - 

_  n-1  n-r-1  n-1  n-r  n-r 

app*^  Var(G),  so  that  we  obtain 


p  ♦  op  r-1 


1,2. 


(44.3) 


This  is  the  form  of  the  autocorrelation  function  of  the  ARMA(1,1)  process. 
Note,  too,  that  the  choice  of  the  slightly  different  structure  for  (44.1)  and 
(4.2)  has  endowed  the  two  parameters  (p,o)  with  physical  significance. 
Choosing  o-O  effectively  sets  B^  to  zero,  so  that  the  process  is  now  simply 
the  BGAR(1)  process,  and  we  can  see  from  (4.3)  that  P)((f’)  becomes  p'',  as 
expected.  If  we  choose  p-0,  then  Y^  becomes  G^  and  {X^l  is  the  BGMA(1) 
process,  and  (4.3)  reduces  to  o/(1*a),  as  it  should.  The  reduction  of  the 
mixed  model  to  its  simpler  forms  when  the  parameters  vanish  is  a  consequence 
of  the  structure  we  have  chosen. 


4.3.  Joint  Distributions  in  the  BGARMA(1,1)  Process. 

The  joint  distribution  of  two  consecutive  observations  of  the  process 
{X  }  can  be  derived  in  the  form  of  the  corresponding  Laplace  transform. 
Thus, 


L(u,v)  -  E{exp(-uX  ^  -vX  ) 
n+ 1  n 

-  E[exp{-u(A  Y  +A*G  )  -  uB  ^  -  vY  ,  -  vB  G  )] 

^  n  n-1  n  n  n^1  n+^l  n-1  n  n 


Lg^(u).E{L^(v+uA)}.E{LQ(uA'  vB) } . 


The  first  two  terms  of  this  product  have  already  been  evaluated  and  we  now 
consider  the  third.  By  considering  appropriate  series  expansions,  we  can 
evaluate  the  third  term  in  the  form 


OD  q» 


y  y  (~u)  .  (-v)  .  r(a0*n)r(pO^m)r(0->-m»n)r(0) 

m-0  n-0  r(o0)r(p0)r(0*n)r(0+m) 

where  0  -  k/(1+a).  Hence,  we  can  show  that  L(u,v)  is  given  by 

-uv 


1 

|0p 

I  1 

(1+u)(1+v) 

[l  ♦u+v 

,F,(0p;Oa,0,y^j 


where  is  the  Hypergeometric  function,  defined  by 


r(a*n)r(b*n)r(c)  z 


... _ _ _n 

2‘Va,b,c,z)  -  r(a)r(b)r(c*n)  n!' 


(4.5) 


The  behavior  and  properties  of  this  function  are  detailed  in  Abramowitz  and 
Stegun  (1964,  Ch.  15).  It  is  easily  verified  that  when  p-0  or  a-0,  the 
appropriate  forms  of  L(u,v)  result  from  (4.5).  The  transform  corresponding 
to  higher-dimensional  distributions  are  more  difficult  to  obtain  in  closed 
form,  although,  in  principle,  series  expansions  of  the  form  of  (4.4)  can  be 


Further,  the  symmetry  of  L(u,v)  in  u  and  v  implies  that  the  conditional 

distributions  of  X  given  X  and  X  given  X  are  identical.  In 

n+1  n  n  n+i 

particular,  we  can  recover  the  conditional  moments  from  (^4.5)  and  it  is 
found  that  regression  is  linear  and  the  conditional  variance  is  quadratic  in 

X  . 

n 

il.'J.  Higher  Order  Mixed  Processes. 

Higher  order  BGARMA  processes  can  be  derived  by  suitable  extensions  of 
the  BGARMA(1,1).  In  particular,  it  is  straightforward  to  construct  a 
BGARMA(1,q)  process  by  replacing  (*4.1)  by  an  MA(q)  form  as  in  (3.6).  Thus, 

q-1 

X  -  I  B  .G  .  +  Y  , 
n  n-i  n-x  n-q' 

replaces  ( *< .  1  )  and  (^<.2)  is  as  before.  The  more  general  problem  of 
extending  to  higher  order  AR  forms  is  more  difficult.  One  way  of  achieving 
it,  however,  is  to  use  mixtures  (random  indexing).  For  details,  the  reader 
is  referred  to  Lewis  (1985). 


5.  DUAL  AND  ITERATED  GAMMA  PROCESSES 


The  first-order  autoregressive  Beta-Gamma  processes  given  in  equation 
(2.1)  has  been  shown  to  be  time-reversible.  This  can  be  a  handicap  in 
modelling  phenomena  such  as  water  run-offs,  which  tend  to  have  'runs  down' 
in  their  sample  paths.  This  is  modelled,  as  noted,  in  a  defective  way  by 
the  GAR(1)  process  of  Gaver  and  Lewis  (1980).  We,  therefore,  look  for  other 
Gamma  processes,  possibly  with  more  than  one  parameter  to  model  dependency 
structure,  which  broaden  the  BGAR(l)  process. 

The  first  process  to  consider  is  the  dual  of  the  BGAR(1).  The  duality 
refers  to  the  fact  that  where  in  (2.1)  we  decrease  the  shape  parameter,  k, 
of  X^_^  by  using  the  Beta-Gamma  transform  and  then  bring  it  up  to  k  by 
adding  an  independent  Gamma  variable,  we  now  increase  k  in  X^_^  by  adding  an 
independent  Gamma  variate  and  then  decrease  the  parameter  to  k  by  using  the 
Beta-Gamma  transform.  Thus,  we  have 


However,  it  can  be  shown  that  the  joint  transform  of  and  is 
(1*u)  ^(1+v)  ^(1+u+v)^  *^2^1  ^  ^  ^  ^ process  is 
time-reversible.  We  thus  have  nothing  new  by  way  of  broadening  the  BGAR(1) 
process . 

Another  approach  to  broadening  the  BGAR(1)  structure  is  to  iterate  the 
process.  Thus,  in  (2.1),  the  left  hand  side  is  a  Gamma(k,B)  random  variable 
and  the  procedure  in  (2.1)  can  be  reapplied.  However,  a  time-reversible 
process  is  again  obtained.  A  better  way  to  iterate  is  to  apply  the  GAR(1) 
procedure  to  (2.1)  and  obtain  a  combination  of  the  BGAR(1)  and  GAR( 1 ) 
processes: 

X  (k,6)  -  Y|B’(kp,kp)X„  ,  ♦  Y^(kp,B)}  ♦  GI  (kY.B)  (5.2) 
n  n  n-i  n  n 

-  YB^(kp,kp)X^_^  ♦  YY^(kp,B)  ♦  GI^(kY,B),  (5.3) 

where  0  i  Y  S  1,  0  ipSI,  p-Y««1,k>0  and  {GI^(kY,6)}  is  a  sequence 
of  i.i.d.  Gamma  innovation  random  variables  with  Laplace-Stieltjes  transform 
(  (  6  ♦Ys  )  /  (  B  +  s ) }  *^ ,  independent  of  {B^}  and  {Y^(kp,B)).  The  condition  that  p 
and  Y  do  not  both  equal  one  is  necessary  to  obtain  an  ergodic  process.  Note 
that  (5.2)  is  different  from  the  combination  given  in  Lawrance  and  Lewis 
(1982)  and  we  denote  it  by  GBGAR(l). 

Now  in  (5.2),  the  case  Y»1  gives  the  BGAR(1)  process,  Y«0  and/or  p-0 
gives  an  i.i.d.  sequence  (X^)  while  p-1  gives  the  GAR(1)  process.  Thus,  we 
should  find  sample  path  behavior  running  from  time-reversibility  to  ’runs 
down'  behavior.  Also,  the  process  is  Markovian  and  has  serial  correlation 

p(r)  -  (Yp)  I’'!  r  -  0,±1  ,+2,  •••.  (5.^0 
For  the  joint  Laplace-Stieltjes  transform  of  X  and  X  , ,  we  have 


Ly  V  (u,v)  -  E{exp(-uX  -vX  ,)} 

X  ,x  .  n  n*i 

n’  n-1 

-  E[exp{-uYB^(kp,kp)X^_^  -  uYY^(kp,B)  -  uGI^(kY,B)  - 


Thus,  it  can  be  seen  that  the  process  is  not  time-reversible  unless  Y-1  (the 

BGAR(l)  case).  The  regression  of  X  on  X  x  is 

n  n-1 

E(Xn|Xn-,-x)  -  p(1)x  +  {1-p(1)}E(X). 

which  is  linear  in  x,  but  it  is  not  the  same  as  E(X  , |X  -x). 

n-1  '  n 

To  separately  Identify  the  parameters  Y  and  p  in  this  process,  one  must 
go  beyond  second-order  properties  of  the  process.  This  is  because  the 
parameters  enter  into  the  correlation  (5.4)  as  a  product.  Higher  order 
residual  analyses  (Lawrance  and  Lewis,  1986)  and  maximum  likelihood 
estimation  will  be  considered  elsewhere. 

6.  NEGATIVE  CORREUTION  AND  NON-LINEAR  PROCESSES 

All  of  the  processes  described  above  are  limited  by  their  serial 
correlations  being  non-negative.  There  are  a  number  of  ways  of  extending 
the  processes  to  give  negative  valued  serial  correlations  and  we  discuss  one 
of  them  in  some  detail.  All  methods  involve  non-linear  functions  of,  say, 
X^_^  in  a  first-order  autoregressive  process.  This  is  necessary  because  of 
the  non-negativity  and  lack  of  symmetry  of  Gamma  disributed  variates. 

6.1.  Antithetic  Variates. 

Let  X  be  a  continuous  random  variable  with  c.d.f.  and  inverse 

c.d.f.  F^'(o),  0  <  a  <1.  Then  the  random  variable  X*  -  {l-Fj^(X)}  is 

called  the  antithetic  variable  to  X.  For  symmetric  two-sided  random 
variables  centered  at  zero,  X*  -  -X.  For  positive  valued  variables  such  as 
Gammas,  X*  has  the  maximum  attainable  negative  correlation  for  bivariate 


Gamma  pairs  (Moran,  1967).  In  particular,  if  k-l  (Exponential), 
X*  -  -Ind-e  *),  but  if  k*1  the  transformation  is  difficult  to  compute. 

If  in  (2.1)  is  replaced  by  X*_^  in  (2.1),  then  a  very  non-linear, 

Markovian  first-order  autoregressive  process  is  obtained.  Serial 
correlations  beyond  lag  one  are  difficult  to  compute. 

6.2.  Coupl ing. 

Caver  and  Lewis  (1980)  introduced  a  scheme  in  the  context  of  the  GAR(1) 
processes  for  cross-coupling  two  Gamma  processes  so  that  the  marginal 
processes  will  have  negative  serial  correlations.  It  is  actually  easier  to 
implement  this  scheme  for  the  BGAR(1)  process  than  for  the  GAR(1)  process, 
because  the  random.  Beta  distributed  coefficients  are  continuous.  We  do  not 
pursue  this  here. 

6.2.  Inverse  Processes. 

A  direct  scheme  for  obtaining  negative  correlation  in  a  Gamma  process  is 
now  given.  It  is  a  generalization  of  a  scheme  given  by  Lewis  (1983)  to 
generate  negatively  correlated  bivariate  Gamma  pairs.  Its  utility  lies  in 
the  fact  that  the  sequence  can  be  generated  with  nothing  but  i.i.d.  Gamma 
variates,  no  numerical  inversions  of  Inverse  distribution  functions  are 
required . 

Thus,  let  B^(k;q-k),  for  q  >  k,  be  a  sequence  of  independent 

Beta(k;q-k)  variates,  independent  of  G^(q)  anjt,  GJ^(k*q ) ,  which  are 

independent  sequences  of  independent  Gamma  variates,  n»l,2,**«.  Also,  let 

X„(k)  be  a  Gamma(k)  variate,  where  k  >  0.  The  idea  is  that  we  want  X  (k)  to 
u  n 

be  small  when  X^_^(k)  is  large,  while  retaining  the  Gamma(k)  marginal 
structure  of  the  process.  We  have 

G;(q) 

X^(k)  -  B^(k,q-k)  - (k)*G'(T)  q  2  k,  n  -  1,2,....  (6.1) 

n-1  n 

Note  that  the  ratio  is  Beta(q;k)  and,  by  the  Beta-Gamma  transform 
(1.7),  the  product  of  this  ratio  with  the  independent  Gamma(k*q)  variable  is 


1 

I 
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a  Gainma(q)  variate.  The  multiplier  B^(k:q-k)  reduces  the  shape  parameter 
from  q  to  k. 

To  obtain  p(1),  the  correlation  between  X  (k)  and  X  ,(k),  we  need 

n  n~i 

I  G'(q)X^_,(k)  1 


k  (kiai  -I  I 

q  B  {X„.,(k)  >  G;(q)l' 


(6.2) 


To  evaluate  the  remaining  expectation  in  (6.2),  we  use  the  fact  that  in  an 

expression  such  as  X  ,(k)/{X  ,(k)  ♦  G'(q)},  the  denominator  is  independent 

n~i  n“i  n 

of  the  ratio.  Then  we  have,  after  some  manipulation,  that 


I  kg 

Xn.i(k)  ♦  G;;(q)l  “  6(k-'-q+1)* 


(6.3) 


Combining  (6.2)  and  (6.3).  we  get  the  suprisingly  simple  result  that 


Corr(Xn.X„.,)  .  ,(t)  -  - 


q  >  k, 


(6.14) 


This  correlation  is  always  negative  and  if  k  >  1  (the  Exponential  case), 
p(1)  has  a  minimum  attainable  value  of  -1/3  when  q  ■  k.  This  is  about 
halfway  to  the  minimum  attainable  correlation  for  bivariate  Exponential 
variables  of  -0.61. 

The  scheme  can  be  Iterated  to  achieve  greater  negative  serial 
correlation.  In  this  and  the  scheme  (6.1),  serial  correlations  of  higher 
order  are  difficult  to  obtain.  However,  since  the  process  is  Markovian,  the 
decay  of  the  absolute  values  of  the  serial  correlations  is  geometrically 
bounded . 


;  -ft' 
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7.  AN  ANALYSIS  OF  TIMES  BETWEEN  FAILURES  OF  A  COMPUTER  SYSTEM 


Hugus  (1982)  used  the  Beta-Gamma  process  to  analyze  a  long  sequence  of 

wind  speeds.  This  sequence  is  very  non-stationary ,  containing  yearly 

cycles.  The  model  actually  used  is  u(n)G*,  where  u(n)  is  a  log-linear 

function  of  n  and  G*  is  a  unit-mean  BGAR(l)  process. 

n 

A  simpler,  stationary  series  of  times  between  failures  of  a  computer 
system  is  analyzed  here.  Although  this  data  is  known  to  be  generated  by  a 
branching  Poisson  process  (Lewis,  196*1),  the  Gamma  model  is  much  simpler  and 
much  more  tractable  than  the  branching  Poisson  process  and  provides  an 
adequate  representation  for  most  purposes.  Modelling  of  these  times  between 
failures  is  important  because,  for  example,  they  represent  times  at  which 
requests  for  service  to  the  computer  are  made. 

In  Figure  7.1,  we  give  a  Gamma  probability  plot  for  the  256  times 
between  failures.  The  fit  appears  adequate,  but  the  goodness  of  fit 
statistics  in  the  table  in  Figure  7.1  must  be  used  with  caution,  since  the 
data  is  serially  correlated  and  two  parameters  have  been  estimated  from  the 
data.  The  parameter  k  (ALPHA  in  the  table)  is  estimated  as  k  ■  0.70**.  Note 
the  two  outliers  in  the  Gamma  probability  plot.  Since  these  are  serially 
adjacent,  they  probably  represent  a  lapse  in  recording  of  computer  failures. 

In  Figure  7.2,  we  show  a  correlation  analysis  of  the  data.  The  decay  in 
correlation  from  ?5(1)  »  0.353  to  higher  lags  is  consistent  with  first-order 
autoregressive  correlation  structure.  The  bands  in  the  figure  are 
approximate  confidence  intervals  for  each  Ji(k)  under  the  assumption  the  true 
correlation  is  zero  for  lag  greater  than  k.  (See  Box  and  Jenkins,  1976,  p. 
35  for  details).  Clearly  the  times  between  failures  are  correlated  and  thus 
a  renewal  model,  say,  for  these  times  between  failures  would  be  inadequate. 

The  partial  autocorrelation  plot  in  Figure  7.2  also  confirms  the  first- 

order  autoregressive  nature  of  the  data.  The  last  panel  in  Figure  7.2  shows 

the  autocorrelation  function  for  the  estimated  residuals,  R  -X  t-p(1)X 

n  n  n-i 

there  is  no  significant  correlation  in  this  series. 

Finally,  in  Figure  7.3.  we  give  empirical  density  functions  (kernel 

density  estimates)  for  the  successive  differences,  ,  n-2,3.’**  and 

successive  ratios  X  /X  ,  n*2,3»»»,  which  were  discussed  in  Section  2.7. 

n  n-1 
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The  differences  show  a  highly  symmetric  and  long-tailed  density  function 
which  is  consistent  with  an  t-Laplace  distribution.  Note  that  the  median  of 
the  differences  is  estimated  to  be  -9.  Given  the  range  of  the  differences, 
this  is  probably  not  significantly  different  from  the  value  of  zero  which 
would  hold  for  the  time-reversible  Beta-Gamma  process.  The  ratios,  » 
have  an  estimated  mean  of  10.812  with  estimated  standard  deviation  of 
^2 . H32 /( 255)  ^  •  2.66.  Thus,  there  is  nothing  here  to  suggest  any 
inadequacy  in  the  Beta-Gamma  model  for  characterizing  the  data. 
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Figure  7. 1 

Gamma  probability  plot  for  times 
between  failures  of  a  computer  system 
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Figure  7.2 

j  Autocorrelation  structure  of  times  between 

I  failures  of  a  computer  system 
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Figure  7.3 


Empirical  density  functions  for  ratios  and 
differences  of  successive  times  between  failures 
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